Genetic analyses

Accompanying ‘AE_TEMPLATE’

General Setup

Create a new analysis directories.

- general directory

- for plots

- for output of summary results

- for baseline tables

- for genetic analyses

- for Cox regression results

source("scripts/functions.R")

source("scripts/pack01.packages.R")

* General packages...

Today = format(as.Date(as.POSIXlt(Sys.time())), "%Y%m%d")
Today.Report = format(as.Date(as.POSIXlt(Sys.time())), "%A, %B %d, %Y")

source("scripts/colors.R")

Background

Collaboration to study PCSK9 in relation to atherosclerotic plaques characteristics.

• Genes.xlsx - list of genes of interest.
• Variants.xlsx - list of variant(s) of interest.

library(openxlsx)

gene_list_df <- read.xlsx(paste0(TARGET_loc, "/targets.xlsx"), sheet = "Genes")

gene_list <- unlist(gene_list_df$Gene)
gene_list

[1] "CXCL10" "PCSK9"  "COL4A1" "COL4A2" "COL3A"  "COL2A"  "LDLR"   "CD36"  

variant_list <- read.xlsx(paste0(TARGET_loc, "/targets.xlsx"), sheet = "Variants")

DT::datatable(variant_list)

NA

This notebook

In this notebook we setup the files for the SNP analysis.

Athero-Express Biobank Study

The Athero-Express Biobank Study (AE) contains plaque material of patients that underwent endarterectomyat two Dutch tertiary referral centers. Details of the study design were described before. Briefly, blood and plaque material were obtained during endarterectomy and stored at -80 ℃. Only carotid endarterectomy (CEA) patients were included in the present study. All patients provided informed consent and the study was approved by the medical ethics committee.

Athero-Express Genomics Study

DNA isolation and genotyping

We genotyped the AE in three separate, but consecutive experiments. In short, DNA was extracted from EDTA blood or (when no blood was available) plaque samples of 1,858 consecutive patients from the Athero-Express Biobank Study and genotyped in 3 batches.

For the Athero-Express Genomics Study 1 (AEGS1) 891 patients (602 males, 262 females, 27 unknown sex), included between 2002 and 2007, were genotyped (440,763 markers) using the Affymetrix Genome-Wide Human SNP Array 5.0 (SNP5) chip (Affymetrix Inc., Santa Clara, CA, USA) at Eurofins Genomics (formerly known as AROS).

For the Athero-Express Genomics Study 2 (AEGS2) 954 patients (640 makes, 313 females, 1 unknown sex), included between 2002 and 2013, were genotyped (587,351 markers) using the Affymetrix AxiomⓇ GW CEU 1 Array (AxM) at the Genome Analysis Center.

For the Athero-Express Genomics Study 3 (AEGS3) 658 patients (448 males, 203 females, 5 unknown sex), included between 2002 and 2016, were genotyped (693,931 markers) using the Illumina GSA MD v1 BeadArray (GSA) at Human Genomics Facility, HUGE-F.

All experiments were carried out according to OECD standards.

Genotyping calling

We used the genotyping calling algorithms as advised by Affymetrix (AEGS1 and AEGS2) and Illumina (AEGS3):

• AEGS1: BRLMM-P
• AEGS2: AxiomGT1
• AEGS3: Illumina GenomeStudio

Quality control after genotyping

After genotype calling, we adhered to community standard quality control and assurance (QCA) procedures of the genotype data from AEGS1, AEGS2, and AEGS3. Samples with low average genotype calling and sex discrepancies (compared to the clinical data available) were excluded. The data was further filtered on:

1. individual (sample) call rate > 97%,
2. SNP call rate > 97%,
3. minor allele frequencies (MAF) > 3%,
4. average heterozygosity rate ± 3.0 s.d.,
5. relatedness (pi-hat > 0.20),
6. Hardy–Weinberg Equilibrium (HWE p < 1.0×10^{−3<>), and}
7. Monomorphic SNPs (< 1.0×10^−6<>).

After QCA 2,493 samples remained, 108 of non-European descent/ancestry, and 156 related pairs. These comprise 890 samples and 407,712 SNPs in AEGS1, 869 samples and 534,508 SNPs in AEGS2, and 649954 samples and 534,508 SNPs in AEGS3 remained.

Imputation

Before phasing using SHAPEIT2, data was lifted to genome build b37 using the liftOver tool from UCSC (https://genome.ucsc.edu/cgi-bin/hgLiftOver). Finally, data was imputed with 1000G phase 3, version 5 and HRC release 1.1 as a reference using the Michigan Imputation Server. These results were further integrated using QCTOOL v2, where HRC imputed variants are given precendence over 1000G phase 3 imputed variants.

Quality control after imputation

We compared quality of the three AEGS datasets, and listed some variables of interest.

• sample type (EDTA blood or plaque)
• genotyping chip used
• reason for filtering

We checked the studytype (AE or not), and identity-by-descent (IBD) within and between datasets to aid in sample mixups, duplicate sample use, and relatedness. In addition, during genotyping quality control samples were identified that deviated from Hardy-Weinberg Equilibrium (HWE), had discordance in sex-coding and genotype sex, and deviated from the principal component analysis (PCA) plot.

We will load the Athero-Express Biobank Study data, and all the samples that were send for genotyping and the final QC’ed sampleList.

Loading data

Loading Athero-Express Biobank Study clinical and biobank data, as well as the SampleList of genetic data. We simply load the previously saved RDS-file and extract the clinical data from that.

cat("* get Athero-Express Biobank Study Database...\n")

* get Athero-Express Biobank Study Database...

AEDB.CEA <- readRDS(file = paste0(OUT_loc, "/20230614.",TRAIT_OF_INTEREST,".AEDB.CEA.RDS"))
AEDB.CEA[1:10, 1:10]

dim(AEDB.CEA)

[1] 2595 1201

AEDB.full <- readRDS(file = paste0(OUT_loc, "/20230614.",TRAIT_OF_INTEREST,".AEDB.FULL.RDS"))
AEDB.full[1:10, 1:10]

dim(AEDB.full)

[1] 3673 1201

cat("\n* get Athero-Express Genomics Study keys...\n")

* get Athero-Express Genomics Study keys...

AEGS123.sampleList.keytable <- fread(paste0(AEGSQC_loc, "/QC/SELECTIONS/20200419.QC.AEGS123.sampleList.keytable.txt"))

dim(AEGS123.sampleList.keytable)

[1] 2124   25

# AEGS123.sampleList.keytable[1:10, 1:10]

Athero-Express Genomics Study

Prepare baseline

Let’s combine the full Athero-Express Biobank Study with the key-table containing the AEGS data.

NOTE: this should sum to 2,124 samples with genotypes.

AEGS <- merge(AEDB.full, AEGS123.sampleList.keytable, by.x = "STUDY_NUMBER", by.y = "STUDY_NUMBER", sort = FALSE,
                  all = TRUE)

dim(AEGS)

[1] 3791 1225

AEGS$UPID.y <- NULL
names(AEGS)[names(AEGS) == "UPID.x"] <- "UPID"
AEGS$Age.y <- NULL
names(AEGS)[names(AEGS) == "Age.x"] <- "Age"

table(AEGS$CHIP, useNA = "ifany")

AffyAxiomCEU     AffySNP5       IllGSA         <NA> 
         918          687          519         1667 

AEGS$GWAS <- AEGS$CHIP
AEGS$GWAS[is.na(AEGS$GWAS)] <- "not genotyped"
AEGS$GWAS[AEGS$GWAS != "not genotyped"] <- "genotyped"

table(AEGS$CHIP, AEGS$GWAS, useNA = "ifany")

              
               genotyped not genotyped
  AffyAxiomCEU       918             0
  AffySNP5           687             0
  IllGSA             519             0
  <NA>                 0          1667

Also a visualisation of the AEGS with AEDB overlaps.

library(UpSetR)
require(ggplot2)
require(plyr)
require(gridExtra)
require(grid)

AEDB.availGWAS = list(
AEGS1 = subset(AEGS, CHIP == "AffySNP5", select = c("STUDY_NUMBER"))[,1],
AEGS2 = subset(AEGS, CHIP == "AffyAxiomCEU", select = c("STUDY_NUMBER"))[,1],
AEGS3 = subset(AEGS, CHIP == "IllGSA", select = c("STUDY_NUMBER"))[,1],
AEDB = AEGS$STUDY_NUMBER)

p1 <- UpSetR::upset(fromList(AEDB.availGWAS), 
                    sets = c("AEDB", "AEGS1", "AEGS2", "AEGS3"), 
                    main.bar.color = c(uithof_color[15], uithof_color[3], uithof_color[2], uithof_color[21]), 
                    mainbar.y.label = "intersection sample size", 
                    sets.bar.color = c(uithof_color[15], uithof_color[2], uithof_color[3], uithof_color[21]), 
                    sets.x.label = "sample size", keep.order = TRUE)

p1

pdf(paste0(PLOT_loc, "/", Today, ".overlap.AEDB_AEGS123.UpSetR.pdf"))
  p1

dev.off()

quartz_off_screen 
                2 

png(paste0(PLOT_loc, "/", Today, ".overlap.AEDB_AEGS123.UpSetR.png"))

  p1
dev.off()

quartz_off_screen 
                2 

rm(p1)

Please refer to the AEDB.CEA.baseline.html document for the details on the informed consent selection.

table(AEGS$Artery_summary, AEGS$QC2018_FILTER)

                                                                                         
                                                                                          family_discard family_keep issue passed
  No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA              0           0     0      0
  carotid (left & right)                                                                              20          21    40   1814
  femoral/iliac (left, right or both sides)                                                            1           0     0     99
  other carotid arteries (common, external)                                                            0           0     0      9
  carotid bypass and injury (left, right or both sides)                                                0           0     0      1
  aneurysmata (carotid & femoral)                                                                      0           0     0      0
  aorta                                                                                                0           0     0      0
  other arteries (renal, popliteal, vertebral)                                                         0           0     0      1
  femoral bypass, angioseal and injury (left, right or both sides)                                     0           0     0      0

table(AEGS$informedconsent, AEGS$QC2018_FILTER)

                                                                                                 
                                                                                                  family_discard family_keep issue passed
  missing                                                                                                      0           0     0      0
  no, died                                                                                                     0           0     0      0
  yes                                                                                                         17          16    24   1411
  yes, health treatment when possible                                                                          2           2    10    312
  yes, no health treatment                                                                                     2           2     2     91
  yes, no health treatment, no commercial business                                                             0           0     0     15
  yes, no tissue, no commerical business                                                                       0           0     0      0
  yes, no tissue, no questionnaires, no medical info, no commercial business                                   0           0     0      0
  yes, no questionnaires, no health treatment, no commercial business                                          0           0     0      1
  yes, no questionnaires, health treatment when possible                                                       0           0     0      2
  yes, no tissue, no questionnaires, no health treatment, no commerical business                               0           0     0      0
  yes, no health treatment, no medical info, no commercial business                                            0           0     3     10
  yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business              0           0     0      0
  yes, no questionnaires, no health treatment                                                                  0           0     0      0
  yes, no tissue, no health treatment                                                                          0           0     0      0
  yes, no tissue, no questionnaires                                                                            0           0     0      0
  yes, no tissue, health treatment when possible                                                               0           0     0      0
  yes, no tissue                                                                                               0           0     0      0
  yes, no commerical business                                                                                  0           1     0     35
  yes, health treatment when possible, no commercial business                                                  0           0     0     25
  yes, no medical info, no commercial business                                                                 0           0     0      4
  yes, no questionnaires                                                                                       0           0     0      1
  yes, no tissue, no questionnaires, no health treatment, no medical info                                      0           0     0      0
  yes, no tissue, no questionnaires, no health treatment, no commercial business                               0           0     0      0
  yes, no medical info                                                                                         0           0     1      5
  yes, no questionnaires, no commercial business                                                               0           0     0      0
  yes, no questionnaires, no health treatment, no medical info                                                 0           0     0      1
  yes, no questionnaires, health treatment when possible, no commercial business                               0           0     0      0
  yes,  no health treatment, no medical info                                                                   0           0     0      5
  no, doesn't want to                                                                                          0           0     0      0
  no, unable to sign                                                                                           0           0     0      0
  no, no reaction                                                                                              0           0     0      0
  no, lost                                                                                                     0           0     0      0
  no, too old                                                                                                  0           0     0      0
  yes, no medical info, health treatment when possible                                                         0           0     0      2
  no (never asked for IC because there was no tissue)                                                          0           0     0      0
  yes, no medical info, no commercial business, health treatment when possible                                 0           0     0      2
  no, endpoint                                                                                                 0           0     0      0
  wil niets invullen, wel alles gebruiken                                                                      0           0     0      0
  second informed concents: yes, no commercial business                                                        0           0     0      2
  nooit geincludeerd                                                                                           0           0     0      0
  yes, not outside EU                                                                                          0           0     0      0
  yes, no DNA                                                                                                  0           0     0      0

AEGSselect <- subset(AEGS, 
                     informedconsent != "missing" & # we are really strict in selecting based on 'informed consent'!
                     informedconsent != "no, died" & 
                     informedconsent != "yes, no tissue, no commerical business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" & 
                     informedconsent != "yes, no tissue, no health treatment" & 
                     informedconsent != "yes, no tissue, no questionnaires" & 
                     informedconsent != "yes, no tissue, health treatment when possible" & 
                     informedconsent != "yes, no tissue" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" & 
                     informedconsent != "no, doesn't want to" & 
                     informedconsent != "no, unable to sign" & 
                     informedconsent != "no, no reaction" & 
                     informedconsent != "no, lost" & 
                     informedconsent != "no, too old" & 
                     informedconsent != "no (never asked for IC because there was no tissue)" & 
                     informedconsent != "no, endpoint" & 
                     informedconsent != "wil niets invullen, wel alles gebruiken" & 
                     informedconsent != "nooit geincludeerd" & 
                     informedconsent != "yes, no DNA")

AEGSselect.CEA <- subset(AEGS, !is.na(QC2018_FILTER) & QC2018_FILTER != "issue" & QC2018_FILTER != "family_discard" &
                     (Artery_summary == "carotid (left & right)" | Artery_summary == "other carotid arteries (common, external)") & # we only want carotids
                     informedconsent != "missing" & # we are really strict in selecting based on 'informed consent'!
                     informedconsent != "no, died" & 
                     informedconsent != "yes, no tissue, no commerical business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no medical info, no commercial business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commerical business" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info, no commercial business" & 
                     informedconsent != "yes, no tissue, no health treatment" & 
                     informedconsent != "yes, no tissue, no questionnaires" & 
                     informedconsent != "yes, no tissue, health treatment when possible" & 
                     informedconsent != "yes, no tissue" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no medical info" & 
                     informedconsent != "yes, no tissue, no questionnaires, no health treatment, no commercial business" & 
                     informedconsent != "no, doesn't want to" & 
                     informedconsent != "no, unable to sign" & 
                     informedconsent != "no, no reaction" & 
                     informedconsent != "no, lost" & 
                     informedconsent != "no, too old" & 
                     informedconsent != "no (never asked for IC because there was no tissue)" & 
                     informedconsent != "no, endpoint" & 
                     informedconsent != "wil niets invullen, wel alles gebruiken" & 
                     informedconsent != "nooit geincludeerd" & 
                     informedconsent != "yes, no DNA")

dim(AEGSselect)

[1] 3673 1224

table(AEGSselect$Artery_summary, AEGSselect$QC2018_FILTER)

                                                                                         
                                                                                          family_discard family_keep issue passed
  No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA              0           0     0      0
  carotid (left & right)                                                                              20          21    40   1814
  femoral/iliac (left, right or both sides)                                                            1           0     0     99
  other carotid arteries (common, external)                                                            0           0     0      9
  carotid bypass and injury (left, right or both sides)                                                0           0     0      1
  aneurysmata (carotid & femoral)                                                                      0           0     0      0
  aorta                                                                                                0           0     0      0
  other arteries (renal, popliteal, vertebral)                                                         0           0     0      1
  femoral bypass, angioseal and injury (left, right or both sides)                                     0           0     0      0

table(AEGSselect$Artery_summary, AEGSselect$CHIP)

                                                                                         
                                                                                          AffyAxiomCEU AffySNP5 IllGSA
  No artery known (yet), no surgery (patient ill, died, exited study), re-numbered to AAA            0        0      0
  carotid (left & right)                                                                           860      550    485
  femoral/iliac (left, right or both sides)                                                          0       72     28
  other carotid arteries (common, external)                                                          2        5      2
  carotid bypass and injury (left, right or both sides)                                              0        1      0
  aneurysmata (carotid & femoral)                                                                    0        0      0
  aorta                                                                                              0        0      0
  other arteries (renal, popliteal, vertebral)                                                       0        1      0
  femoral bypass, angioseal and injury (left, right or both sides)                                   0        0      0

table(AEGSselect$QC2018_FILTER, AEGSselect$CHIP)

                
                 AffyAxiomCEU AffySNP5 IllGSA
  family_discard           12        6      3
  family_keep               8        1     12
  issue                    37        1      2
  passed                  805      621    498

table(AEGSselect$QC2018_FILTER, AEGSselect$SAMPLE_TYPE)

                
                 EDTA_blood plaque unknown
  family_discard         15      6       0
  family_keep            13      8       0
  issue                  26     14       0
  passed               1200    723       1

AEDB.temp <- subset(AEGSselect,  select = c("STUDY_NUMBER", "UPID", "Age", "Gender", "Hospital", "Artery_summary", "QC2018_FILTER", "CHIP", "SAMPLE_TYPE"))
require(labelled)
AEDB.temp$Gender <- to_factor(AEDB.temp$Gender)
AEDB.temp$Hospital <- to_factor(AEDB.temp$Hospital)
AEDB.temp$Artery_summary <- to_factor(AEDB.temp$Artery_summary)
AEDB.temp$QC2018_FILTER <- to_factor(AEDB.temp$QC2018_FILTER)
AEDB.temp$CHIP <- to_factor(AEDB.temp$CHIP)
AEDB.temp$SAMPLE_TYPE <- to_factor(AEDB.temp$SAMPLE_TYPE)

DT::datatable(AEDB.temp[1:10,], caption = "Excerpt of the whole AEDB.", rownames = FALSE)

rm(AEDB.temp)

Athero-Express Genomics Study Baseline Characteristics

Showing the baseline table of the Athero-Express Genomics Study.

```r
# Create baseline tables
# http://rstudio-pubs-static.s3.amazonaws.com/13321_da314633db924dc78986a850813a50d5.html
AEGSselect$GWAS <- to_factor(AEGSselect$GWAS)
AEGSselect$CHIP <- to_factor(AEGSselect$CHIP)
AEGSselect$PCA <- to_factor(AEGSselect$PCA)
AEGSselect$SAMPLE_TYPE <- to_factor(AEGSselect$SAMPLE_TYPE)
AEGSselect$informedconsent <- to_factor(AEGSselect$informedconsent)

AEGSselect.CEA$GWAS <- to_factor(AEGSselect.CEA$GWAS)
AEGSselect.CEA$CHIP <- to_factor(AEGSselect.CEA$CHIP)
AEGSselect.CEA$PCA <- to_factor(AEGSselect.CEA$PCA)
AEGSselect.CEA$SAMPLE_TYPE <- to_factor(AEGSselect.CEA$SAMPLE_TYPE)
AEGSselect.CEA$informedconsent <- to_factor(AEGSselect.CEA$informedconsent)


cat(\===========================================================================================\n\)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxuIn0= -->

===========================================================================================

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuY2F0KFxcQ1JFQVRFIEJBU0VMSU5FIFRBQkxFXFxuXFwpXG5gYGBcbmBgYCJ9 -->

```r
```r
cat(\CREATE BASELINE TABLE\n\)

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiQ1JFQVRFIEJBU0VMSU5FIFRBQkxFXG4ifQ== -->

CREATE BASELINE TABLE

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuIyBCYXNlbGluZSB0YWJsZSB2YXJpYWJsZXNcbmJhc2V0YWJsZV92YXJzID0gYyhcXEhvc3BpdGFsXFwsIFxuICAgICAgICAgICAgICAgICAgIFxcQWdlXFwsIFxcR2VuZGVyXFwsIFxuICAgICAgICAgICAgICAgICAgIFxcVENfZmluYWxcXCwgXFxMRExfZmluYWxcXCwgXFxIRExfZmluYWxcXCwgXFxUR19maW5hbFxcLCBcbiAgICAgICAgICAgICAgICAgICBcXHN5c3RvbGljXFwsIFxcZGlhc3RvbGlcXCwgXFxHRlJfTURSRFxcLCBcXEJNSVxcLCBcbiAgICAgICAgICAgICAgICAgICBcXEtET1FJXFwsIFxcQk1JX1dIT1xcLCBcbiAgICAgICAgICAgICAgICAgICBcXFNtb2tlckN1cnJlbnRcXCwgXFxlQ2lnYXJldHRlc1xcLCBcXGVQYWNrWWVhcnNTbW9raW5nXFwsXG4gICAgICAgICAgICAgICAgICAgXFxEaWFiZXRlc1N0YXR1c1xcLCBcXEh5cGVydGVuc2lvbi5zZWxmcmVwb3J0XFwsIFxcSHlwZXJ0ZW5zaW9uLnNlbGZyZXBvcnRkcnVnXFwsIFxcSHlwZXJ0ZW5zaW9uLmNvbXBvc2l0ZVxcLCBcbiAgICAgICAgICAgICAgICAgICBcXEh5cGVydGVuc2lvbi5kcnVnc1xcLCBcXE1lZC5hbnRpY29hZ3VsYW50c1xcLCBcXE1lZC5hbGwuYW50aXBsYXRlbGV0XFwsIFxcTWVkLlN0YXRpbi5MTERcXCwgXG4gICAgICAgICAgICAgICAgICAgXFxTdHJva2VfRHhcXCwgXFxzeW1wdFxcLCBcXFN5bXB0b21zLjVHXFwsIFxccmVzdGVub3NcXCxcbiAgICAgICAgICAgICAgICAgICBcXEVQX2NvbXBvc2l0ZVxcLCBcXEVQX2NvbXBvc2l0ZV90aW1lXFwsXG4gICAgICAgICAgICAgICAgICAgXFxtYWNtZWFuMFxcLCBcXHNtY21lYW4wXFwsIFxcTWFjcm9waGFnZXMuYmluXFwsIFxcU01DLmJpblxcLCBcXG5ldXRyb3BoaWxzXFwsIFxcTWFzdF9jZWxsc19wbGFxdWVcXCwgXFx2ZXNzZWxfZGVuc2l0eV9hdmVyYWdlZFxcLFxuICAgICAgICAgICAgICAgICAgIFxcSVBILmJpblxcLCBcbiAgICAgICAgICAgICAgICAgICBcXENhbGMuYmluXFwsIFxcQ29sbGFnZW4uYmluXFwsIFxuICAgICAgICAgICAgICAgICAgIFxcRmF0LmJpbl8xMFxcLCBcXEZhdC5iaW5fNDBcXCwgXFxPdmVyYWxsUGxhcXVlUGhlbm90eXBlXFwsIFxcUGxhcXVlX1Z1bG5lcmFiaWxpdHlfSW5kZXhcXCxcbiAgICAgICAgICAgICAgICAgICBcXFNNQ19yYW5rTm9ybVxcLCBcXE1BQ19yYW5rTm9ybVxcLCBcXE5ldXRyb3BoaWxzX3JhbmtOb3JtXFwsIFxcTWFzdENlbGxzX3JhbmtOb3JtXFwsIFxcVmVzc2VsRGVuc2l0eV9yYW5rTm9ybVxcLFxuICAgICAgICAgICAgICAgICAgIFxcR1dBU1xcLCBcXENISVBcXCwgXFxQQ0FcXCxcbiAgICAgICAgICAgICAgICAgICBcXFBDU0s5X3BsYXNtYVxcLCBcXFBDU0s5X3BsYXNtYV9yYW5rTm9ybVxcKVxuXG5iYXNldGFibGVfYmluID0gYyhcXEdlbmRlclxcLCBcbiAgICAgICAgICAgICAgICAgIFxcS0RPUUlcXCwgXFxCTUlfV0hPXFwsIFxuICAgICAgICAgICAgICAgICAgXFxTbW9rZXJDdXJyZW50XFwsIFxuICAgICAgICAgICAgICAgICAgXFxEaWFiZXRlc1N0YXR1c1xcLCBcXEh5cGVydGVuc2lvbi5zZWxmcmVwb3J0XFwsIFxcSHlwZXJ0ZW5zaW9uLnNlbGZyZXBvcnRkcnVnXFwsIFxcSHlwZXJ0ZW5zaW9uLmNvbXBvc2l0ZVxcLCBcbiAgICAgICAgICAgICAgICAgIFxcSHlwZXJ0ZW5zaW9uLmRydWdzXFwsIFxcTWVkLmFudGljb2FndWxhbnRzXFwsIFxcTWVkLmFsbC5hbnRpcGxhdGVsZXRcXCwgXFxNZWQuU3RhdGluLkxMRFxcLCBcbiAgICAgICAgICAgICAgICAgIFxcU3Ryb2tlX0R4XFwsIFxcc3ltcHRcXCwgXFxTeW1wdG9tcy41R1xcLCBcXHJlc3Rlbm9zXFwsXG4gICAgICAgICAgICAgICAgICBcXEVQX2NvbXBvc2l0ZVxcLCBcXE1hY3JvcGhhZ2VzLmJpblxcLCBcXFNNQy5iaW5cXCxcbiAgICAgICAgICAgICAgICAgIFxcSVBILmJpblxcLCBcbiAgICAgICAgICAgICAgICAgIFxcQ2FsYy5iaW5cXCwgXFxDb2xsYWdlbi5iaW5cXCwgXG4gICAgICAgICAgICAgICAgICBcXEZhdC5iaW5fMTBcXCwgXFxGYXQuYmluXzQwXFwsIFxcT3ZlcmFsbFBsYXF1ZVBoZW5vdHlwZVxcLCBcXFBsYXF1ZV9WdWxuZXJhYmlsaXR5X0luZGV4XFwsXG4gICAgICAgICAgICAgICAgICBcXEdXQVNcXCwgXFxDSElQXFwsIFxcUENBXFwpXG5cbmJhc2V0YWJsZV9iaW5cbmBgYFxuYGBgIn0= -->

```r
```r
# Baseline table variables
basetable_vars = c(\Hospital\, 
                   \Age\, \Gender\, 
                   \TC_final\, \LDL_final\, \HDL_final\, \TG_final\, 
                   \systolic\, \diastoli\, \GFR_MDRD\, \BMI\, 
                   \KDOQI\, \BMI_WHO\, 
                   \SmokerCurrent\, \eCigarettes\, \ePackYearsSmoking\,
                   \DiabetesStatus\, \Hypertension.selfreport\, \Hypertension.selfreportdrug\, \Hypertension.composite\, 
                   \Hypertension.drugs\, \Med.anticoagulants\, \Med.all.antiplatelet\, \Med.Statin.LLD\, 
                   \Stroke_Dx\, \sympt\, \Symptoms.5G\, \restenos\,
                   \EP_composite\, \EP_composite_time\,
                   \macmean0\, \smcmean0\, \Macrophages.bin\, \SMC.bin\, \neutrophils\, \Mast_cells_plaque\, \vessel_density_averaged\,
                   \IPH.bin\, 
                   \Calc.bin\, \Collagen.bin\, 
                   \Fat.bin_10\, \Fat.bin_40\, \OverallPlaquePhenotype\, \Plaque_Vulnerability_Index\,
                   \SMC_rankNorm\, \MAC_rankNorm\, \Neutrophils_rankNorm\, \MastCells_rankNorm\, \VesselDensity_rankNorm\,
                   \GWAS\, \CHIP\, \PCA\,
                   \PCSK9_plasma\, \PCSK9_plasma_rankNorm\)

basetable_bin = c(\Gender\, 
                  \KDOQI\, \BMI_WHO\, 
                  \SmokerCurrent\, 
                  \DiabetesStatus\, \Hypertension.selfreport\, \Hypertension.selfreportdrug\, \Hypertension.composite\, 
                  \Hypertension.drugs\, \Med.anticoagulants\, \Med.all.antiplatelet\, \Med.Statin.LLD\, 
                  \Stroke_Dx\, \sympt\, \Symptoms.5G\, \restenos\,
                  \EP_composite\, \Macrophages.bin\, \SMC.bin\,
                  \IPH.bin\, 
                  \Calc.bin\, \Collagen.bin\, 
                  \Fat.bin_10\, \Fat.bin_40\, \OverallPlaquePhenotype\, \Plaque_Vulnerability_Index\,
                  \GWAS\, \CHIP\, \PCA\)

basetable_bin

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiIFsxXSBcXEdlbmRlclxcICAgICAgICAgICAgICAgICAgICAgIFxcS0RPUUlcXCAgICAgICAgICAgICAgICAgICAgICAgXFxCTUlfV0hPXFwgICAgICAgICAgICAgICAgICAgIFxuIFs0XSBcXFNtb2tlckN1cnJlbnRcXCAgICAgICAgICAgICAgIFxcRGlhYmV0ZXNTdGF0dXNcXCAgICAgICAgICAgICAgXFxIeXBlcnRlbnNpb24uc2VsZnJlcG9ydFxcICAgIFxuIFs3XSBcXEh5cGVydGVuc2lvbi5zZWxmcmVwb3J0ZHJ1Z1xcIFxcSHlwZXJ0ZW5zaW9uLmNvbXBvc2l0ZVxcICAgICAgXFxIeXBlcnRlbnNpb24uZHJ1Z3NcXCAgICAgICAgIFxuWzEwXSBcXE1lZC5hbnRpY29hZ3VsYW50c1xcICAgICAgICAgIFxcTWVkLmFsbC5hbnRpcGxhdGVsZXRcXCAgICAgICAgXFxNZWQuU3RhdGluLkxMRFxcICAgICAgICAgICAgIFxuWzEzXSBcXFN0cm9rZV9EeFxcICAgICAgICAgICAgICAgICAgIFxcc3ltcHRcXCAgICAgICAgICAgICAgICAgICAgICAgXFxTeW1wdG9tcy41R1xcICAgICAgICAgICAgICAgIFxuWzE2XSBcXHJlc3Rlbm9zXFwgICAgICAgICAgICAgICAgICAgIFxcRVBfY29tcG9zaXRlXFwgICAgICAgICAgICAgICAgXFxNYWNyb3BoYWdlcy5iaW5cXCAgICAgICAgICAgIFxuWzE5XSBcXFNNQy5iaW5cXCAgICAgICAgICAgICAgICAgICAgIFxcSVBILmJpblxcICAgICAgICAgICAgICAgICAgICAgXFxDYWxjLmJpblxcICAgICAgICAgICAgICAgICAgIFxuWzIyXSBcXENvbGxhZ2VuLmJpblxcICAgICAgICAgICAgICAgIFxcRmF0LmJpbl8xMFxcICAgICAgICAgICAgICAgICAgXFxGYXQuYmluXzQwXFwgICAgICAgICAgICAgICAgIFxuWzI1XSBcXE92ZXJhbGxQbGFxdWVQaGVub3R5cGVcXCAgICAgIFxcUGxhcXVlX1Z1bG5lcmFiaWxpdHlfSW5kZXhcXCAgXFxHV0FTXFwgICAgICAgICAgICAgICAgICAgICAgIFxuWzI4XSBcXENISVBcXCAgICAgICAgICAgICAgICAgICAgICAgIFxcUENBXFwgICAgICAgICAgICAgICAgICAgICAgICBcbiJ9 -->

[1]     _WHO 
[4]     .selfreport 
[7] .selfreportdrug .composite  .drugs 
[10] .anticoagulants  .all.antiplatelet  .Statin.LLD 
[13] _Dx    .5G 
[16]   _composite  .bin 
[19] .bin  .bin  .bin 
[22] .bin  .bin_10  .bin_40 
[25]   _Vulnerability_Index   
[28]    

<!-- rnb-output-end -->

<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuYGBgclxuYmFzZXRhYmxlX2NvbiA9IGJhc2V0YWJsZV92YXJzWyFiYXNldGFibGVfdmFycyAlaW4lIGJhc2V0YWJsZV9iaW5dXG5iYXNldGFibGVfY29uXG5gYGBcbmBgYCJ9 -->

```r
```r
basetable_con = basetable_vars[!basetable_vars %in% basetable_bin]
basetable_con

<!-- rnb-source-end -->

<!-- rnb-output-begin eyJkYXRhIjoiIFsxXSBcXEhvc3BpdGFsXFwgICAgICAgICAgICAgICAgXFxBZ2VcXCAgICAgICAgICAgICAgICAgICAgIFxcVENfZmluYWxcXCAgICAgICAgICAgICAgICBcXExETF9maW5hbFxcICAgICAgICAgICAgICBcbiBbNV0gXFxIRExfZmluYWxcXCAgICAgICAgICAgICAgIFxcVEdfZmluYWxcXCAgICAgICAgICAgICAgICBcXHN5c3RvbGljXFwgICAgICAgICAgICAgICAgXFxkaWFzdG9saVxcICAgICAgICAgICAgICAgXG4gWzldIFxcR0ZSX01EUkRcXCAgICAgICAgICAgICAgICBcXEJNSVxcICAgICAgICAgICAgICAgICAgICAgXFxlQ2lnYXJldHRlc1xcICAgICAgICAgICAgIFxcZVBhY2tZZWFyc1Ntb2tpbmdcXCAgICAgIFxuWzEzXSBcXEVQX2NvbXBvc2l0ZV90aW1lXFwgICAgICAgXFxtYWNtZWFuMFxcICAgICAgICAgICAgICAgIFxcc21jbWVhbjBcXCAgICAgICAgICAgICAgICBcXG5ldXRyb3BoaWxzXFwgICAgICAgICAgICBcblsxN10gXFxNYXN0X2NlbGxzX3BsYXF1ZVxcICAgICAgIFxcdmVzc2VsX2RlbnNpdHlfYXZlcmFnZWRcXCBcXFNNQ19yYW5rTm9ybVxcICAgICAgICAgICAgXFxNQUNfcmFua05vcm1cXCAgICAgICAgICAgXG5bMjFdIFxcTmV1dHJvcGhpbHNfcmFua05vcm1cXCAgICBcXE1hc3RDZWxsc19yYW5rTm9ybVxcICAgICAgXFxWZXNzZWxEZW5zaXR5X3JhbmtOb3JtXFwgIFxcUENTSzlfcGxhc21hXFwgICAgICAgICAgIFxuWzI1XSBcXFBDU0s5X3BsYXNtYV9yYW5rTm9ybVxcICBcbiJ9 -->

[1]     _final  _final 
[5] _final  _final     
[9] _MDRD       
[13] _composite_time       
[17] _cells_plaque  _density_averaged _rankNorm  _rankNorm 
[21] _rankNorm  _rankNorm  _rankNorm  _plasma 
[25] _plasma_rankNorm 

<!-- rnb-output-end -->

<!-- rnb-chunk-end -->


<!-- rnb-text-begin -->


All Athero-Express Genomics Study data (n = 2,011), compared to the *remaining*, \_un_genotyped Athero-Express Biobank
Study.


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuY2F0KFwiXFxuPT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PT09PVxcblwiKVxuY2F0KFwiRElTUExBWSBCQVNFTElORSBUQUJMRVxcblwiKVxuXG5BRUdTc2VsZWN0LnRhYmxlT25lID0gcHJpbnQoQ3JlYXRlVGFibGVPbmUodmFycyA9IGJhc2V0YWJsZV92YXJzLCBcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIyBmYWN0b3JWYXJzID0gYmFzZXRhYmxlX2JpbixcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgc3RyYXRhID0gXCJHV0FTXCIsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIGRhdGEgPSBBRUdTc2VsZWN0LCBpbmNsdWRlTkEgPSBUUlVFKSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgIG5vbm5vcm1hbCA9IGMoKSwgbWlzc2luZyA9IFRSVUUsXG4gICAgICAgICAgICAgICAgICAgICAgICAgIHF1b3RlID0gRkFMU0UsIG5vU3BhY2VzID0gRkFMU0UsIHNob3dBbGxMZXZlbHMgPSBUUlVFLCBleHBsYWluID0gVFJVRSwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgIGZvcm1hdCA9IFwicGZcIiwgXG4gICAgICAgICAgICAgICAgICAgICAgICAgIGNvbnREaWdpdHMgPSAzKVssMTo2XVxuYGBgIn0= -->

```r
cat("\n===========================================================================================\n")
cat("DISPLAY BASELINE TABLE\n")

AEGSselect.tableOne = print(CreateTableOne(vars = basetable_vars, 
                                         # factorVars = basetable_bin,
                                         strata = "GWAS",
                                         data = AEGSselect, includeNA = TRUE), 
                          nonnormal = c(), missing = TRUE,
                          quote = FALSE, noSpaces = FALSE, showAllLevels = TRUE, explain = TRUE, 
                          format = "pf", 
                          contDigits = 3)[,1:6]

Baseline of the valid, CEA and genotyped data.

AEGSselect.CEA.tableOne = print(CreateTableOne(vars = basetable_vars, 
                                         # factorVars = basetable_bin,
                                         strata = "Gender",
                                         data = AEGSselect.CEA, includeNA = TRUE), 
                          nonnormal = c(), missing = TRUE,
                          quote = FALSE, noSpaces = FALSE, showAllLevels = TRUE, explain = TRUE, 
                          format = "pf", 
                          contDigits = 3)[,1:6]

Baseline writing

Let’s save the baseline characteristics of the Athero-Express Genomics Study.

# Write basetable
require(openxlsx)

write.xlsx(file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AEGS.BaselineTable.xlsx"), 
           format(as.data.frame(AEGSselect.tableOne), digits = 5, scientific = FALSE),
           rowNames = TRUE, 
           colNames = TRUE, 
           sheetName = "AEGS_Base_AEDB", overwrite = TRUE)

write.xlsx(file = paste0(BASELINE_loc, "/",Today,".",PROJECTNAME,".AEGS.CEA.BaselineTable.xlsx"), 
           format(as.data.frame(AEGSselect.CEA.tableOne), digits = 5, scientific = FALSE),
           rowNames = TRUE, 
           colNames = TRUE, 
           sheetName = "AEGS_Base_CEA_sex", overwrite = TRUE)

SampleLists

Autosomal data.

We are ready to make a sampleList for use with the imputed data.

require(openxlsx)

temp <- subset(AEGS,
               GWAS == "genotyped",
               select = c("ID_1", "ID_2", "UPID", "STUDY_NUMBER", # ID_2 is the order of samples!
                          "QC2018_FINAL", "QC2018_FILTER", "OriginalOrder_postMichImp_QC",
                          "AEGS_type", "CHIP", "STUDY_TYPE", "SAMPLE_TYPE", "PCA",
                          "PC1", "PC2", "PC3", "PC4", "PC5",
                          "PC6", "PC7", "PC8", "PC9", "PC10",
                          "Sex", "Age", "ORyear", 
                          "Calc.bin", "Collagen.bin", 
                          "Fat.bin_10", "Fat.bin_40", "IPH.bin", 
                          "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm",
                          "Plaque_Vulnerability_Index",
                          "PCSK9_plasma", "PCSK9_plasma_rankNorm")) # Select some phenotype of interest
dim(temp)

[1] 2124   38

# Fix things
attach(temp)

temp[,"Calcification"] <- NA
temp$Calcification[Calc.bin == "no/minor"] <- "control"
temp$Calcification[Calc.bin == "moderate/heavy"] <- "case"

temp[,"Collagen"] <- NA
temp$Collagen[Collagen.bin == "no/minor"] <- "control"
temp$Collagen[Collagen.bin == "moderate/heavy"] <- "case"

temp[,"Fat10"] <- NA
temp$Fat10[Fat.bin_10 == "<10%"] <- "control"
temp$Fat10[Fat.bin_10 == ">10%"] <- "case"

temp[,"Fat40"] <- NA
temp$Fat40[Fat.bin_40 == "<40%"] <- "control"
temp$Fat40[Fat.bin_40 == ">40%"] <- "case"

temp[,"IPH"] <- NA
temp$IPH[IPH.bin == "no"] <- "control"
temp$IPH[IPH.bin == "yes"] <- "case"

temp[,"PVI"] <- NA
temp$PVI[Plaque_Vulnerability_Index == "0"] <- "PVI_cat0"
temp$PVI[Plaque_Vulnerability_Index == "1"] <- "PVI_cat1"
temp$PVI[Plaque_Vulnerability_Index == "2"] <- "PVI_cat2"
temp$PVI[Plaque_Vulnerability_Index == "3"] <- "PVI_cat3"
temp$PVI[Plaque_Vulnerability_Index == "4"] <- "PVI_cat4"
temp$PVI[Plaque_Vulnerability_Index == "5"] <- "PVI_cat5"

temp$Plaque_Vulnerability_Index <- temp$PVI
temp$PVI <- NULL

detach(temp)

# Making selection variable
attach(temp)
temp[,"SELECTION"] <- "not_selected"
temp$SELECTION[(QC2018_FILTER=="passed" | QC2018_FILTER=="family_keep") & (STUDY_TYPE=="CEA" & PCA=="EUR")] <- "selected"
detach(temp)
table(temp$SELECTION, temp$QC2018_FILTER)

              
               family_discard family_keep issue passed
  not_selected             23           0    41    152
  selected                  0          21     0   1887

table(temp$SELECTION, temp$STUDY_TYPE)

              
                CEA  FEA Other
  not_selected   94  109    12
  selected     1908    0     0

table(temp$SELECTION, temp$PCA)

              
                EUR nonEUR
  not_selected  167     45
  selected     1908      0

table(temp$Plaque_Vulnerability_Index)

PVI_cat0 PVI_cat1 PVI_cat2 PVI_cat3 PVI_cat4 PVI_cat5 
     453      348      422      478      221       84 

# Check https://www.well.ox.ac.uk/~gav/snptest/#multinomial_tests for details on multiple categories testing in SNPTEST

AEGS123_sample.list <- temp[order(temp$OriginalOrder_postMichImp_QC),]

AEGS123_sample.list$missing <- 0

sample_file_aegs <- dplyr::select(AEGS123_sample.list,
                                  ID_1, ID_2, missing, # ID_2 is the order of samples - that way we always know what the order should be
                                  UPID, STUDY_NUMBER, 
                                  QC2018_FINAL, QC2018_FILTER, SELECTION,
                                  AEGS_type, CHIP, STUDY_TYPE, SAMPLE_TYPE, PCA,
                                  PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9, PC10,
                                  Sex, Age, ORyear, 
                                  Calcification, Collagen, 
                                  Fat10, Fat40, IPH, 
                                  SMC_rankNorm, MAC_rankNorm, Neutrophils_rankNorm, MastCells_rankNorm, VesselDensity_rankNorm,
                                  Plaque_Vulnerability_Index,
                                  PCSK9_plasma, PCSK9_plasma_rankNorm)  %>%
  mutate_if(is.numeric, as.character) %>%
  mutate(SAMPLE_TYPE = gsub(' ', '_', SAMPLE_TYPE)) %>%
  add_row(.before = 1, 
          ID_1 = "0", ID_2 = "0", missing = "0", 
          UPID = "D", STUDY_NUMBER = "C",
          QC2018_FINAL = "D", QC2018_FILTER = "D", SELECTION = "D",
          AEGS_type = "D", CHIP = "D", STUDY_TYPE = "D", SAMPLE_TYPE = "D", PCA = "D",
          PC1 = "C", PC2 = "C", PC3 = "C", PC4 = "C", PC5 = "C", PC6 = "C", PC7 = "C", PC8 = "C", PC9 = "C", PC10 = "C",
          Sex = "D", Age = "C", ORyear = "C", 
          Calcification = "B", Collagen = "B", 
          Fat10 = "B", Fat40 = "B", IPH = "B", 
          SMC_rankNorm = "P", MAC_rankNorm = "P", Neutrophils_rankNorm = "P", MastCells_rankNorm = "P", VesselDensity_rankNorm = "P",
          Plaque_Vulnerability_Index = "D",
          PCSK9_plasma = "P", PCSK9_plasma_rankNorm = "P") %>% ## identifiers: index for these is 1, and all base variables have 0 as identifier
  print()

dim(sample_file_aegs)

[1] 2125   39

fwrite(sample_file_aegs,
       file = paste0(SNP_loc, "/",Today,".",PROJECTNAME,".AEGS123.sample"),
       na = "NA", sep = "\t", quote = FALSE,
       row.names = FALSE, col.names = TRUE,
       showProgress = TRUE, verbose = TRUE)

This installation of data.table has not been compiled with OpenMP support.
  omp_get_num_procs()            1
  R_DATATABLE_NUM_PROCS_PERCENT  unset (default 50)
  R_DATATABLE_NUM_THREADS        unset
  R_DATATABLE_THROTTLE           unset (default 1024)
  omp_get_thread_limit()         1
  omp_get_max_threads()          1
  OMP_THREAD_LIMIT               unset
  OMP_NUM_THREADS                1
  RestoreAfterFork               true
  data.table is using 1 threads with throttle==1024. See ?setDTthreads.
No list columns are present. Setting sep2='' otherwise quote='auto' would quote fields containing sep2.
Column writers: 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 
args.doRowNames=0 args.rowNames=0 doQuote=0 args.nrow=2125 args.ncol=39 eolLen=1
maxLineLen=924. Found in 0.000s
Writing bom (false), yaml (0 characters) and column names (true) ... done in 0.001s
Writing 2125 rows in 1 batches of 2125 rows (each buffer size 8MB, showProgress=1, nth=1)

require(DT)
DT::datatable(sample_file_aegs, caption = "AEGS: final sample list of genotyped AE patients after quality control.", rownames = FALSE)

Warning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlWarning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

rm(temp, temp2, temp3)

Warning: object 'temp2' not foundWarning: object 'temp3' not found

Females only

This is the selection for females only.

require(openxlsx)

temp <- subset(AEGS,
               GWAS == "genotyped",
               select = c("ID_1", "ID_2", "UPID", "STUDY_NUMBER", # ID_2 is the order of samples!
                          "QC2018_FINAL", "QC2018_FILTER", "OriginalOrder_postMichImp_QC",
                          "AEGS_type", "CHIP", "STUDY_TYPE", "SAMPLE_TYPE", "PCA",
                          "PC1", "PC2", "PC3", "PC4", "PC5",
                          "PC6", "PC7", "PC8", "PC9", "PC10",
                          "Sex", "Age", "ORyear", 
                          "BMI",
                          "Calc.bin", "Collagen.bin", 
                          "Fat.bin_10", "Fat.bin_40", "IPH.bin", 
                          "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm", 
                          "OverallPlaquePhenotype", "Plaque_Vulnerability_Index",
                          "PCSK9_plasma", "PCSK9_plasma_rankNorm")) # Select some phenotype of interest
dim(temp)

[1] 2124   40

# Fix things
attach(temp)

temp[,"Calcification"] <- NA
temp$Calcification[Calc.bin == "no/minor"] <- "control"
temp$Calcification[Calc.bin == "moderate/heavy"] <- "case"

temp[,"Collagen"] <- NA
temp$Collagen[Collagen.bin == "no/minor"] <- "control"
temp$Collagen[Collagen.bin == "moderate/heavy"] <- "case"

temp[,"Fat10"] <- NA
temp$Fat10[Fat.bin_10 == "<10%"] <- "control"
temp$Fat10[Fat.bin_10 == ">10%"] <- "case"

temp[,"Fat40"] <- NA
temp$Fat40[Fat.bin_40 == "<40%"] <- "control"
temp$Fat40[Fat.bin_40 == ">40%"] <- "case"

temp[,"IPH"] <- NA
temp$IPH[IPH.bin == "no"] <- "control"
temp$IPH[IPH.bin == "yes"] <- "case"

temp[,"PVI"] <- NA
temp$PVI[Plaque_Vulnerability_Index == "0"] <- "PVI_cat0"
temp$PVI[Plaque_Vulnerability_Index == "1"] <- "PVI_cat1"
temp$PVI[Plaque_Vulnerability_Index == "2"] <- "PVI_cat2"
temp$PVI[Plaque_Vulnerability_Index == "3"] <- "PVI_cat3"
temp$PVI[Plaque_Vulnerability_Index == "4"] <- "PVI_cat4"
temp$PVI[Plaque_Vulnerability_Index == "5"] <- "PVI_cat5"

temp$Plaque_Vulnerability_Index <- temp$PVI
temp$PVI <- NULL

detach(temp)

# Making selection variable
attach(temp)
temp[,"SELECTION"] <- "not_selected"
temp$SELECTION[(QC2018_FILTER=="passed" | QC2018_FILTER=="family_keep") & (STUDY_TYPE=="CEA" & PCA=="EUR") & Sex=="F"] <- "selected"
detach(temp)
table(temp$SELECTION, temp$QC2018_FILTER)

              
               family_discard family_keep issue passed
  not_selected             23          15    41   1427
  selected                  0           6     0    612

table(temp$SELECTION, temp$STUDY_TYPE)

              
                CEA  FEA Other
  not_selected 1384  109    12
  selected      618    0     0

table(temp$SELECTION, temp$PCA)

              
                EUR nonEUR
  not_selected 1457     45
  selected      618      0

AEGS123_sample.list <- temp[order(temp$OriginalOrder_postMichImp_QC),]

AEGS123_sample.list$missing <- 0

sample_file_aegsF <- dplyr::select(AEGS123_sample.list,
                                  ID_1, ID_2, missing, # ID_2 is the order of samples - that way we always know what the order should be
                                  UPID, STUDY_NUMBER, 
                                  QC2018_FINAL, QC2018_FILTER, SELECTION,
                                  AEGS_type, CHIP, STUDY_TYPE, SAMPLE_TYPE, PCA,
                                  PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9, PC10,
                                  Sex, Age, ORyear, 
                                  Calcification, Collagen, 
                                  Fat10, Fat40, IPH, 
                                  SMC_rankNorm, MAC_rankNorm, Neutrophils_rankNorm, MastCells_rankNorm, VesselDensity_rankNorm,
                                  Plaque_Vulnerability_Index,
                                  PCSK9_plasma, PCSK9_plasma_rankNorm)  %>%
  mutate_if(is.numeric, as.character) %>%
  mutate(SAMPLE_TYPE = gsub(' ', '_', SAMPLE_TYPE)) %>%
  add_row(.before = 1, 
          ID_1 = "0", ID_2 = "0", missing = "0", 
          UPID = "D", STUDY_NUMBER = "C",
          QC2018_FINAL = "D", QC2018_FILTER = "D", SELECTION = "D",
          AEGS_type = "D", CHIP = "D", STUDY_TYPE = "D", SAMPLE_TYPE = "D", PCA = "D",
          PC1 = "C", PC2 = "C", PC3 = "C", PC4 = "C", PC5 = "C", PC6 = "C", PC7 = "C", PC8 = "C", PC9 = "C", PC10 = "C",
          Sex = "D", Age = "C", ORyear = "C", 
          Calcification = "B", Collagen = "B", 
          Fat10 = "B", Fat40 = "B", IPH = "B", 
          SMC_rankNorm = "P", MAC_rankNorm = "P", Neutrophils_rankNorm = "P", MastCells_rankNorm = "P", VesselDensity_rankNorm = "P",
          Plaque_Vulnerability_Index = "D",
          PCSK9_plasma = "P", PCSK9_plasma_rankNorm = "P") %>% ## identifiers: index for these is 1, and all base variables have 0 as identifier
  print()

dim(sample_file_aegsF)

[1] 2125   39

fwrite(sample_file_aegsF,
       file = paste0(SNP_loc, "/",Today,".",PROJECTNAME,".AEGS123.females.sample"),
       na = "NA", sep = "\t", quote = FALSE,
       row.names = FALSE, col.names = TRUE,
       showProgress = TRUE, verbose = TRUE)

This installation of data.table has not been compiled with OpenMP support.
  omp_get_num_procs()            1
  R_DATATABLE_NUM_PROCS_PERCENT  unset (default 50)
  R_DATATABLE_NUM_THREADS        unset
  R_DATATABLE_THROTTLE           unset (default 1024)
  omp_get_thread_limit()         1
  omp_get_max_threads()          1
  OMP_THREAD_LIMIT               unset
  OMP_NUM_THREADS                1
  RestoreAfterFork               true
  data.table is using 1 threads with throttle==1024. See ?setDTthreads.
No list columns are present. Setting sep2='' otherwise quote='auto' would quote fields containing sep2.
Column writers: 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 
args.doRowNames=0 args.rowNames=0 doQuote=0 args.nrow=2125 args.ncol=39 eolLen=1
maxLineLen=924. Found in 0.000s
Writing bom (false), yaml (0 characters) and column names (true) ... done in 0.001s
Writing 2125 rows in 1 batches of 2125 rows (each buffer size 8MB, showProgress=1, nth=1)

require(DT)
DT::datatable(sample_file_aegsF, caption = "AEGS: final sample list of genotyped AE patients after quality control.", rownames = FALSE)

Warning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlWarning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

rm(temp)

Males only

This is the selection for males only.

require(openxlsx)

temp <- subset(AEGS,
               GWAS == "genotyped",
               select = c("ID_1", "ID_2", "UPID", "STUDY_NUMBER", # ID_2 is the order of samples!
                          "QC2018_FINAL", "QC2018_FILTER", "OriginalOrder_postMichImp_QC",
                          "AEGS_type", "CHIP", "STUDY_TYPE", "SAMPLE_TYPE", "PCA",
                          "PC1", "PC2", "PC3", "PC4", "PC5",
                          "PC6", "PC7", "PC8", "PC9", "PC10",
                          "Sex", "Age", "ORyear", 
                          "BMI",
                          "Calc.bin", "Collagen.bin", 
                          "Fat.bin_10", "Fat.bin_40", "IPH.bin", 
                          "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm", 
                          "OverallPlaquePhenotype", "Plaque_Vulnerability_Index",
                          "PCSK9_plasma", "PCSK9_plasma_rankNorm")) # Select some phenotype of interest
dim(temp)

[1] 2124   40

# Fix things
attach(temp)

temp[,"Calcification"] <- NA
temp$Calcification[Calc.bin == "no/minor"] <- "control"
temp$Calcification[Calc.bin == "moderate/heavy"] <- "case"

temp[,"Collagen"] <- NA
temp$Collagen[Collagen.bin == "no/minor"] <- "control"
temp$Collagen[Collagen.bin == "moderate/heavy"] <- "case"

temp[,"Fat10"] <- NA
temp$Fat10[Fat.bin_10 == "<10%"] <- "control"
temp$Fat10[Fat.bin_10 == ">10%"] <- "case"

temp[,"Fat40"] <- NA
temp$Fat40[Fat.bin_40 == "<40%"] <- "control"
temp$Fat40[Fat.bin_40 == ">40%"] <- "case"

temp[,"IPH"] <- NA
temp$IPH[IPH.bin == "no"] <- "control"
temp$IPH[IPH.bin == "yes"] <- "case"

temp[,"PVI"] <- NA
temp$PVI[Plaque_Vulnerability_Index == "0"] <- "PVI_cat0"
temp$PVI[Plaque_Vulnerability_Index == "1"] <- "PVI_cat1"
temp$PVI[Plaque_Vulnerability_Index == "2"] <- "PVI_cat2"
temp$PVI[Plaque_Vulnerability_Index == "3"] <- "PVI_cat3"
temp$PVI[Plaque_Vulnerability_Index == "4"] <- "PVI_cat4"
temp$PVI[Plaque_Vulnerability_Index == "5"] <- "PVI_cat5"

temp$Plaque_Vulnerability_Index <- temp$PVI
temp$PVI <- NULL

detach(temp)

# Making selection variable
attach(temp)
temp[,"SELECTION"] <- "not_selected"
temp$SELECTION[(QC2018_FILTER=="passed" | QC2018_FILTER=="family_keep") & (STUDY_TYPE=="CEA" & PCA=="EUR") & Sex=="M"] <- "selected"
detach(temp)
table(temp$SELECTION, temp$QC2018_FILTER)

              
               family_discard family_keep issue passed
  not_selected             23           6    41    764
  selected                  0          15     0   1275

table(temp$SELECTION, temp$STUDY_TYPE)

              
                CEA  FEA Other
  not_selected  712  109    12
  selected     1290    0     0

table(temp$SELECTION, temp$PCA)

              
                EUR nonEUR
  not_selected  785     45
  selected     1290      0

AEGS123_sample.list <- temp[order(temp$OriginalOrder_postMichImp_QC),]

AEGS123_sample.list$missing <- 0

sample_file_aegsM <- dplyr::select(AEGS123_sample.list,
                                  ID_1, ID_2, missing, # ID_2 is the order of samples - that way we always know what the order should be
                                  UPID, STUDY_NUMBER, 
                                  QC2018_FINAL, QC2018_FILTER, SELECTION,
                                  AEGS_type, CHIP, STUDY_TYPE, SAMPLE_TYPE, PCA,
                                  PC1, PC2, PC3, PC4, PC5, PC6, PC7, PC8, PC9, PC10,
                                  Sex, Age, ORyear, 
                                  Calcification, Collagen, 
                                  Fat10, Fat40, IPH, 
                                  SMC_rankNorm, MAC_rankNorm, Neutrophils_rankNorm, MastCells_rankNorm, VesselDensity_rankNorm,
                                  Plaque_Vulnerability_Index,
                                  PCSK9_plasma, PCSK9_plasma_rankNorm)  %>%
  mutate_if(is.numeric, as.character) %>%
  mutate(SAMPLE_TYPE = gsub(' ', '_', SAMPLE_TYPE)) %>%
  add_row(.before = 1, 
          ID_1 = "0", ID_2 = "0", missing = "0", 
          UPID = "D", STUDY_NUMBER = "C",
          QC2018_FINAL = "D", QC2018_FILTER = "D", SELECTION = "D",
          AEGS_type = "D", CHIP = "D", STUDY_TYPE = "D", SAMPLE_TYPE = "D", PCA = "D",
          PC1 = "C", PC2 = "C", PC3 = "C", PC4 = "C", PC5 = "C", PC6 = "C", PC7 = "C", PC8 = "C", PC9 = "C", PC10 = "C",
          Sex = "D", Age = "C", ORyear = "C", 
          Calcification = "B", Collagen = "B", 
          Fat10 = "B", Fat40 = "B", IPH = "B", 
          SMC_rankNorm = "P", MAC_rankNorm = "P", Neutrophils_rankNorm = "P", MastCells_rankNorm = "P", VesselDensity_rankNorm = "P",
          Plaque_Vulnerability_Index = "D",
          PCSK9_plasma = "P", PCSK9_plasma_rankNorm = "P") %>% ## identifiers: index for these is 1, and all base variables have 0 as identifier
  print()

dim(sample_file_aegsM)

[1] 2125   39

fwrite(sample_file_aegsM,
       file = paste0(SNP_loc, "/",Today,".",PROJECTNAME,".AEGS123.males.sample"),
       na = "NA", sep = "\t", quote = FALSE,
       row.names = FALSE, col.names = TRUE,
       showProgress = TRUE, verbose = TRUE)

This installation of data.table has not been compiled with OpenMP support.
  omp_get_num_procs()            1
  R_DATATABLE_NUM_PROCS_PERCENT  unset (default 50)
  R_DATATABLE_NUM_THREADS        unset
  R_DATATABLE_THROTTLE           unset (default 1024)
  omp_get_thread_limit()         1
  omp_get_max_threads()          1
  OMP_THREAD_LIMIT               unset
  OMP_NUM_THREADS                1
  RestoreAfterFork               true
  data.table is using 1 threads with throttle==1024. See ?setDTthreads.
No list columns are present. Setting sep2='' otherwise quote='auto' would quote fields containing sep2.
Column writers: 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 
args.doRowNames=0 args.rowNames=0 doQuote=0 args.nrow=2125 args.ncol=39 eolLen=1
maxLineLen=924. Found in 0.000s
Writing bom (false), yaml (0 characters) and column names (true) ... done in 0.001s
Writing 2125 rows in 1 batches of 2125 rows (each buffer size 8MB, showProgress=1, nth=1)

require(DT)
DT::datatable(sample_file_aegsM, caption = "AEGS: final sample list of genotyped AE patients after quality control.", rownames = FALSE)

Warning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlWarning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

rm(temp)

X-chromosome data

The X-chromosome data is taken from previously imputed data based on 1000G phase 3 (version 5) and GoNL5. For some reason, imputing on the Michigan Imputation Server was not successful (ACTION point).

Here we load in the sample files for the three datasets of the X chromosomal data. We should:

• filter out samples that did not pass quality control, ending up with 2,124 sample
• re-order the data to fit the other autosomal data.

AEGS123_chrX <- fread(paste0(MICHIMP_loc, "/_chr23_1kg_gonl5/aegs.raw.1kg_gonl5.chr23.mappings.txt"))
names(AEGS123_chrX)[names(AEGS123_chrX) == "ID_1"] <- "SampleID_postImpChrX"

AEGS123_AllChr <- merge(AEGS123_chrX, sample_file_aegs, by.x = "SampleID_postMichImp", by.y = "ID_1", 
                        all.x = TRUE, 
                        sort = FALSE)

names(AEGS123_AllChr)[names(AEGS123_AllChr) == "ID_2.y"] <- "ID_2"
names(AEGS123_AllChr)[names(AEGS123_AllChr) == "STUDY_TYPE.y"] <- "STUDY_TYPE"
names(AEGS123_AllChr)[names(AEGS123_AllChr) == "SampleID_postMichImp"] <- "ID_1"
AEGS123_AllChr$missing.x <- NULL
AEGS123_AllChr$missing.y <- NULL
AEGS123_AllChr$STUDY_TYPE.x <- NULL
AEGS123_AllChr$ID_2.x <- NULL

dim(AEGS123_AllChr)

[1] 2176   44

str(AEGS123_AllChr)

Classes ‘data.table’ and 'data.frame':  2176 obs. of  44 variables:
 $ ID_1                      : chr  "0" "UPID00126-UPID00126-A318-0034" "UPID01799-UPID01799-A318-0492" "UPID01890-UPID01890-A318-0151" ...
 $ SampleID_postImpChrX      : chr  "0" "UPID00126-UPID00126-A318-0034" "UPID01799-UPID01799-A318-0492" "UPID01890-UPID01890-A318-0151" ...
 $ MappingID                 : chr  "0" "UPID00126-UPID00126-A318-0034" "UPID01799-UPID01799-A318-0492" "UPID01890-UPID01890-A318-0151" ...
 $ FID_forQC                 : chr  "0" "UPID00126" "UPID01799" "UPID01890" ...
 $ IID_forQC                 : chr  "0" "UPID00126-UPID00126-A318-0034" "UPID01799-UPID01799-A318-0492" "UPID01890-UPID01890-A318-0151" ...
 $ SampleID_postQC           : chr  "0" "UPID00126-UPID00126-A318-0034" "UPID01799-UPID01799-A318-0492" "UPID01890-UPID01890-A318-0151" ...
 $ ChrX_Order                : int  0 1 2 NA 3 4 5 6 NA 7 ...
 $ ID_2                      : chr  "0" "22" "325" NA ...
 $ UPID                      : chr  "D" NA NA NA ...
 $ STUDY_NUMBER              : chr  "C" "901" "1164" NA ...
 $ QC2018_FINAL              : chr  "D" NA NA NA ...
 $ QC2018_FILTER             : chr  "D" "passed" "passed" NA ...
 $ SELECTION                 : chr  "D" "selected" "selected" NA ...
 $ AEGS_type                 : chr  "D" "AEGS1" "AEGS1" NA ...
 $ CHIP                      : chr  "D" "AffySNP5" "AffySNP5" NA ...
 $ STUDY_TYPE                : chr  "D" "CEA" "CEA" NA ...
 $ SAMPLE_TYPE               : chr  "D" "plaque" "plaque" NA ...
 $ PCA                       : chr  "D" "EUR" "EUR" NA ...
 $ PC1                       : chr  "C" "-0.0316" "-0.0131" NA ...
 $ PC2                       : chr  "C" "0.0161" "-0.0069" NA ...
 $ PC3                       : chr  "C" "-0.01" "0.0383" NA ...
 $ PC4                       : chr  "C" "-0.0131" "-0.005" NA ...
 $ PC5                       : chr  "C" "0.0059" "0.0078" NA ...
 $ PC6                       : chr  "C" "-0.0283" "0.0243" NA ...
 $ PC7                       : chr  "C" "-2e-04" "-0.018" NA ...
 $ PC8                       : chr  "C" "0.0033" "0.0022" NA ...
 $ PC9                       : chr  "C" "-0.0317" "-0.0154" NA ...
 $ PC10                      : chr  "C" "-0.0054" "0.007" NA ...
 $ Sex                       : chr  "D" "F" "M" NA ...
 $ Age                       : chr  "C" NA NA NA ...
 $ ORyear                    : chr  "C" NA NA NA ...
 $ Calcification             : chr  "B" NA NA NA ...
 $ Collagen                  : chr  "B" NA NA NA ...
 $ Fat10                     : chr  "B" NA NA NA ...
 $ Fat40                     : chr  "B" NA NA NA ...
 $ IPH                       : chr  "B" NA NA NA ...
 $ SMC_rankNorm              : chr  "P" NA NA NA ...
 $ MAC_rankNorm              : chr  "P" NA NA NA ...
 $ Neutrophils_rankNorm      : chr  "P" NA NA NA ...
 $ MastCells_rankNorm        : chr  "P" NA NA NA ...
 $ VesselDensity_rankNorm    : chr  "P" NA NA NA ...
 $ Plaque_Vulnerability_Index: chr  "D" NA NA NA ...
 $ PCSK9_plasma              : chr  "P" NA NA NA ...
 $ PCSK9_plasma_rankNorm     : chr  "P" NA NA NA ...
 - attr(*, ".internal.selfref")=<externalptr> 

This seems fine, let’s filter; we can use this file to filter the genetic data. And we create another file to re-order the data.

AEGS123_AllChrQC <- subset(AEGS123_AllChr,
               !is.na(QC2018_FILTER),
               select = c("ID_1", "ID_2", "UPID", "STUDY_NUMBER", "SampleID_postImpChrX",
                          "QC2018_FINAL", "QC2018_FILTER", "SELECTION",
                          "AEGS_type", "CHIP", "STUDY_TYPE", "SAMPLE_TYPE",
                          "PC1", "PC2", "PC3", "PC4", "PC5",
                          "PC6", "PC7", "PC8", "PC9", "PC10",
                          "Sex", "Age", "ORyear", 
                          "Calcification", "Collagen", 
                          "Fat10", "Fat40", "IPH", 
                          "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm",
                          "Plaque_Vulnerability_Index",
                          "PCSK9_plasma", "PCSK9_plasma_rankNorm"))

AEGS123_AllChrQC_reorder <-AEGS123_AllChrQC[order(AEGS123_AllChrQC$ID_2),] # remember: ID_2 is the order of samples

AEGS123_AllChrQC_filtered <- subset(AEGS123_AllChr,
               is.na(QC2018_FILTER),
               select = c("ID_1", "ID_2", "UPID", "STUDY_NUMBER", "SampleID_postImpChrX",
                          "QC2018_FINAL", "QC2018_FILTER", "SELECTION",
                          "AEGS_type", "CHIP", "STUDY_TYPE", "SAMPLE_TYPE",
                          "PC1", "PC2", "PC3", "PC4", "PC5",
                          "PC6", "PC7", "PC8", "PC9", "PC10",
                          "Sex", "Age", "ORyear", 
                          "Calcification", "Collagen", 
                          "Fat10", "Fat40", "IPH", 
                          "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm",
                          "Plaque_Vulnerability_Index",
                          "PCSK9_plasma", "PCSK9_plasma_rankNorm"))

fwrite(AEGS123_AllChrQC_reorder,
       file = paste0(SNP_loc, "/",Today,".",PROJECTNAME,".AEGS123.chrX.sample"),
       na = "NA", sep = "\t", quote = FALSE,
       row.names = FALSE, col.names = TRUE,
       showProgress = TRUE, verbose = TRUE)

This installation of data.table has not been compiled with OpenMP support.
  omp_get_num_procs()            1
  R_DATATABLE_NUM_PROCS_PERCENT  unset (default 50)
  R_DATATABLE_NUM_THREADS        unset
  R_DATATABLE_THROTTLE           unset (default 1024)
  omp_get_thread_limit()         1
  omp_get_max_threads()          1
  OMP_THREAD_LIMIT               unset
  OMP_NUM_THREADS                1
  RestoreAfterFork               true
  data.table is using 1 threads with throttle==1024. See ?setDTthreads.
No list columns are present. Setting sep2='' otherwise quote='auto' would quote fields containing sep2.
Column writers: 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 12 
args.doRowNames=0 args.rowNames=0 doQuote=0 args.nrow=2013 args.ncol=38 eolLen=1
maxLineLen=991. Found in 0.000s
Writing bom (false), yaml (0 characters) and column names (true) ... done in 0.001s
Writing 2013 rows in 1 batches of 2013 rows (each buffer size 8MB, showProgress=1, nth=1)

require(DT)
DT::datatable(AEGS123_AllChrQC, caption = "AEGS: final sample list of genotyped AE patients after quality control (chromosome X).", rownames = FALSE)

Warning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.htmlWarning: It seems your data is too big for client-side DataTables. You may consider server-side processing: https://rstudio.github.io/DT/server.html

GWASToolKit preparation

VariantLists: top SNPs

Here we create a variantlist.txt file used by GWASToolKit for analysis.These are the variants with LD r² => 0.6 around each of the 15 top SNPs associated with CAC.

variant_list

temp <- subset(variant_list, select = c("VariantID", "Chr", "BP"))

fwrite(temp,
       file = paste0(SNP_loc, "/variantlist.txt"),
       na = "NA", sep = "\t", quote = FALSE,
       row.names = FALSE, col.names = FALSE,
       showProgress = TRUE, verbose = TRUE)

This installation of data.table has not been compiled with OpenMP support.
  omp_get_num_procs()            1
  R_DATATABLE_NUM_PROCS_PERCENT  unset (default 50)
  R_DATATABLE_NUM_THREADS        unset
  R_DATATABLE_THROTTLE           unset (default 1024)
  omp_get_thread_limit()         1
  omp_get_max_threads()          1
  OMP_THREAD_LIMIT               unset
  OMP_NUM_THREADS                1
  RestoreAfterFork               true
  data.table is using 1 threads with throttle==1024. See ?setDTthreads.
No list columns are present. Setting sep2='' otherwise quote='auto' would quote fields containing sep2.
Column writers: 12 5 5 
args.doRowNames=0 args.rowNames=0 doQuote=0 args.nrow=5 args.ncol=3 eolLen=1
maxLineLen=140. Found in 0.000s
Writing bom (false), yaml (0 characters) and column names (false) ... done in 0.002s
Writing 5 rows in 1 batches of 5 rows (each buffer size 8MB, showProgress=1, nth=1)

rm(temp)

Covariates

Here we create a covariates.txt file used by GWASToolKit for analysis.

library(tidyverse)
# for 'overall' analyses
c("Age Sex PC1 PC2 ORyear") %>% write_lines(paste0(SNP_loc, "/covariates.txt"))

# for sex-specific analyses
c("Age PC1 PC2 ORyear") %>% write_lines(paste0(SNP_loc, "/covariates.sex.txt"))

Phenotypes

Here we create a phenotypes.txt file used by GWASToolKit for analysis.

library(tidyverse)
c("Calcification", "Collagen", "Fat10", "Fat40", "IPH", "SMC_rankNorm", "MAC_rankNorm", "Neutrophils_rankNorm", "MastCells_rankNorm", "VesselDensity_rankNorm") %>% write_lines(paste0(SNP_loc, "/phenotypes.txt"))

c("Plaque_Vulnerability_Index") %>% write_lines(paste0(SNP_loc, "/phenotypes.pvi.txt"))

Session information

---

Version:      v1.2.1
Last update:  2023-06-14
Written by:   Sander W. van der Laan (s.w.vanderlaan-2[at]umcutrecht.nl).
Description:  Script to get some Athero-Express Biobank Study baseline characteristics.
Minimum requirements: R version 3.4.3 (2017-06-30) -- 'Single Candle', Mac OS X El Capitan

Changes log
* v1.2.1 Fixed issue with baseline writing.
* v1.2.0 Update to the study database. 
* v1.1.0 Major update to WORCS system. 
* v1.0.6 Small bug fixes.
* v1.0.5 Added png for overlap-figure.
* v1.0.5 Removed obsolete references to objects.
* v1.0.4 Fixed a mistake in the chr X sample-file creation. Now the order matches the chr X data.
* v1.0.3 Fixed weight of files (limit of 10Mb per file for templates). Renamed entire repo.
* v1.0.2 Added sex-specific .sample-files. Added GWASToolKit input-files.
* v1.0.0 Initial version. Add 'plaque vulnerability index', Fixed baseline table, added codes, and results. Created sample-files.

---

sessionInfo()

R version 4.3.0 (2023-04-21)
Platform: x86_64-apple-darwin22.4.0 (64-bit)
Running under: macOS 14.0

Matrix products: default
BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
LAPACK: /usr/local/Cellar/r/4.3.0_1/lib/R/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Amsterdam
tzcode source: internal

attached base packages:
[1] grid      tools     stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gridExtra_2.3        plyr_1.8.8           rmarkdown_2.22       patchwork_1.1.2.9000 labelled_2.11.0      sjPlot_2.8.14        UpSetR_1.4.0        
 [8] ggpubr_0.6.0         forestplot_3.1.1     abind_1.4-5          checkmate_2.2.0      pheatmap_1.0.12      devtools_2.4.5       usethis_2.2.0       
[15] BlandAltmanLeh_0.3.1 tableone_0.13.2      openxlsx_4.2.5.2     haven_2.5.2          eeptools_1.2.5       DT_0.28              knitr_1.43          
[22] lubridate_1.9.2      forcats_1.0.0        stringr_1.5.0        purrr_1.0.1          tibble_3.2.1         ggplot2_3.4.2        tidyverse_2.0.0     
[29] data.table_1.14.8    naniar_1.0.0         tidyr_1.3.0          dplyr_1.1.2          optparse_1.7.3       readr_2.1.4          pander_0.6.5        
[36] R.utils_2.12.2       R.oo_1.25.0          R.methodsS3_1.8.2    worcs_0.1.10         credentials_1.3.2   

loaded via a namespace (and not attached):
  [1] fs_1.6.2               matrixStats_1.0.0      spatstat.sparse_3.0-1  httr_1.4.6             RColorBrewer_1.1-3     insight_0.19.2        
  [7] gh_1.4.0               profvis_0.3.8          sctransform_0.3.5      backports_1.4.1        sjlabelled_1.2.0       utf8_1.2.3            
 [13] R6_2.5.1               lazyeval_0.2.2         uwot_0.1.14            prereg_0.6.0           urlchecker_1.0.1       withr_2.5.0           
 [19] sp_1.6-1               prettyunits_1.1.1      progressr_0.13.0       cli_3.6.1              performance_0.10.4     spatstat.explore_3.2-1
 [25] sandwich_3.0-2         labeling_0.4.2         sass_0.4.6             Seurat_4.3.0           mvtnorm_1.2-1          spatstat.data_3.0-1   
 [31] proxy_0.4-27           ggridges_0.5.4         pbapply_1.7-0          askpass_1.1            sessioninfo_1.2.2      parallelly_1.36.0     
 [37] rstudioapi_0.14        generics_0.1.3         vroom_1.6.3            crosstalk_1.2.0        ica_1.0-3              spatstat.random_3.1-5 
 [43] car_3.1-2              zip_2.3.0              Matrix_1.5-4.1         fansi_1.0.4            lifecycle_1.0.3        multcomp_1.4-23       
 [49] yaml_2.3.7             carData_3.0-5          Rtsne_0.16             promises_1.2.0.1       crayon_1.5.2           miniUI_0.1.1.1        
 [55] lattice_0.21-8         cowplot_1.1.1          sys_3.4.2              pillar_1.9.0           boot_1.3-28.1          estimability_1.4.1    
 [61] future.apply_1.11.0    codetools_0.2-19       leiden_0.4.3           glue_1.6.2             remotes_2.4.2          vcd_1.4-11            
 [67] vctrs_0.6.2            png_0.1-8              gtable_0.3.3           cachem_1.0.8           xfun_0.39              mime_0.12             
 [73] survey_4.2-1           coda_0.19-4            survival_3.5-5         tinytex_0.45           ellipsis_0.3.2         fitdistrplus_1.1-11   
 [79] TH.data_1.1-2          ROCR_1.0-11            nlme_3.1-162           bit64_4.0.5            RcppAnnoy_0.0.20       bslib_0.4.2           
 [85] irlba_2.3.5.1          KernSmooth_2.23-21     colorspace_2.1-0       DBI_1.1.3              tidyselect_1.2.0       processx_3.8.1        
 [91] emmeans_1.8.6          bit_4.0.5              compiler_4.3.0         curl_5.0.0             plotly_4.10.2          bayestestR_0.13.1     
 [97] scales_1.2.1           lmtest_0.9-40          callr_3.7.3            digest_0.6.31          goftest_1.2-3          spatstat.utils_3.0-3  
[103] minqa_1.2.5            htmltools_0.5.5        pkgconfig_2.0.3        rticles_0.25           lme4_1.1-33            fastmap_1.1.1         
[109] rlang_1.1.1            htmlwidgets_1.6.2      shiny_1.7.4            farver_2.1.1           jquerylib_0.1.4        zoo_1.8-12            
[115] jsonlite_1.8.5         magrittr_2.0.3         munsell_0.5.0          Rcpp_1.0.10            reticulate_1.29        visdat_0.6.0          
[121] stringi_1.7.12         MASS_7.3-60            pkgbuild_1.4.0         parallel_4.3.0         listenv_0.9.0          ggrepel_0.9.3         
[127] sjmisc_2.8.9           deldir_1.0-9           ggeffects_1.2.2        splines_4.3.0          tensor_1.5             hms_1.1.3             
[133] sjstats_0.18.2         ps_1.7.5               igraph_1.4.3           ranger_0.15.1          spatstat.geom_3.2-1    ggsignif_0.6.4        
[139] reshape2_1.4.4         pkgload_1.3.2          evaluate_0.21          SeuratObject_4.1.3     mitools_2.4            modelr_0.1.11         
[145] renv_0.17.3            nloptr_2.0.3           tzdb_0.4.0             httpuv_1.6.11          RANN_2.6.1             openssl_2.0.6         
[151] getopt_1.20.3          polyclip_1.10-4        future_1.32.0          scattermore_1.1        broom_1.0.4            xtable_1.8-4          
[157] e1071_1.7-13           rstatix_0.7.2          later_1.3.1            class_7.3-22           viridisLite_0.4.2      gert_1.9.2            
[163] arm_1.13-1             memoise_2.0.1          cluster_2.1.4          timechange_0.2.0       globals_0.16.2        

Saving environment

save.image(paste0(PROJECT_loc, "/",Today,".",PROJECTNAME,".SNP_analyses.RData"))

© 1979-2023 Sander W. van der Laan | s.w.vanderlaan[at]gmail.com | vanderlaan.science. |